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Abstract 

Introduction: Tlie aim of tliis study was to evaluate, for the first time, the differences in gene expression profiles of 
normal and osteoarthritic (OA) subchondral bone in human subjects. 

Methods: Following histological assessment of the integrity of overlying cartilage and the severity of bone 
abnormality by micro-computed tomography, we isolated total RNA from regions of interest from human OA 
(n = 20) and non-OA (n = 5) knee lateral tibial (LT) and medial tibial (MT) plateaus. A whole-genome profiling study 
was performed on an Agilent microarray platform and analyzed using Agilent GeneSpring GXl 1.5. Confirmatory 
quantitative reverse-transcription polymerase chain reaction (qRT-PCR) analysis was performed on samples from 
9 OA individuals to confirm differential expression of 85 genes identified by microarray. Ingenuity Pathway Analysis 
(IPA) was used to investigate canonical pathways and immunohistochemical staining was performed to validate 
protein expression levels in samples. 

Results: A total of 972 differentially expressed genes were identified (fold change > ± 2, P <0.05) between 
LT (minimal degeneration) and MT (significant degeneration) regions from OA samples; these data implicated 
279 canonical pathways in IPA. The qRT-PCR data strongly confirmed the accuracy of microarray results (R^ = 0.58, 
P <0.0001). Novel pathways were identified in this study including Periostin (POSTN) and Leptin (LEP), which are 
implicated in bone remodeling by osteoblasts. 

Conclusions: To the best of our knowledge, this study represents the most comprehensive direct assessment to 
date of gene expression profiling in OA subchondral bone. This study provides insights that could contribute to the 
development of new biomarkers and therapeutic strategies for OA. 



Introduction 

Osteoarthritis (OA) is the most prevalent musculoskeletal 
disorder and the most common form of arthritis among 
older individuals in most countries worldwide; OA consti- 
tutes a large economic burden due to the associated costs 
of medical care and lost wages [1,2]. Although degene- 
ration of cartilage is the major characteristic of OA, the 
disease also involves the entire joint organ, including 
structural modifications of underlying subchondral bone, 
pathological changes of the meniscus and synovitis [3,4]. 
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Maintenance of normal joint structure and function relies 
on load adaptation of the cartilage and bone. Disruption 
of the physiological relationship between these tissues 
contributes to the development of OA pathology [5]. 

Increasing evidence indicates that the subchondral 
bone, consisting of the subchondral plate and subchon- 
dral spongiosa, plays a major role in the initiation and 
progression of OA [6]. Magnetic resonance imaging- 
based visualization of the whole knee structure demon- 
strates that increased tibial subchondral bone volume is 
associated with severity of knee OA [7]. Kraus and col- 
leagues demonstrated that subchondral bone texture can 
be used as a biomarker to predict progression of knee 
OA [8]. A study comparing two guinea pig strains indi- 
cates that an increasing rate of knee subchondral bone 
remodeling is associated with the process of cartilage 
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deterioration [9]. Moreover, human bone cell culture 
studies have shown that factors released from bone cells 
might affect cartilage metabolism [10,11]. These studies 
provide insights into a temporal relationship between 
subchondral structural changes and alterations in articu- 
lar cartilage during the development of OA; they also 
underscore the importance of elucidating the molecular 
changes in human subchondral bone to improve our un- 
derstanding of the pathogenesis of OA. 

Whole-genome microarray is a common technology 
for studying the behavior of many genes simultaneously. 
All of the gene expression microarray profiling studies 
in OA so far have been performed on human articular 
cartilage, meniscus or synovium; however, none have 
been performed on human subchondral bone tissue dir- 
ectly. Only one study has reported the gene expression 
profiles of subchondral bone in an early OA experimen- 
tal mouse model [12]. A few studies have evaluated the 
gene expression profiles of distal trabecular bone from 
human OA [13,14], but this site is remote from sub- 
chondral bone rather than a reflection of alterations 
locally in the OA joint. This paucity of subchondral bone 
microarray studies is most probably due to the difficul- 
ties associated with isolation of high-quality RNA from 
subchondral bone. As described in our previous study 
[15], we have developed a method of precisely harvesting 
specific regions of subchondral bone tissue and for 
subsequently isolating high-quality RNA from these 
specimens. Our methodology has made it possible to 
perform microarray analyses of human subchondral 
bone samples. 

Our goal in this study was to evaluate the association 
of subchondral bone gene expression with bone histo- 
morphometry at sites of intact articular cartilage and 
osteoarthritic lesioned cartilage. To our knowledge, this 
is the first study to successfully perform microarray ana- 
lyses of human knee subchondral bone in OA and non- 
OA samples, thereby providing clues to the pathogenic 
mechanisms of OA that could inform development of 
new diagnostic markers and therapeutic targets. 

Methods 

Human knee joint tissues 

Human osteoarthritic tibial plateaus with medial compart- 
ment dominant OA and macroscopically normal lateral 
compartments were obtained during total knee joint re- 
placement surgery from knee OA patients. The diagnosis 
of OA was based on the criteria for knee OA of the Ameri- 
can College of Rheumatology [16]. Normal tibial plateaus 
were obtained at autopsy or within 8 hours after amputa- 
tion surgery from donors with nonarthritic knee joints. 

To ensure consistency of sampling of prespecified re- 
gions of interest, the anatomic orientation was indicated 
on the freshly isolated specimens by marker pen and 



then all specimens were stored immediately in liquid 
nitrogen. Methods for precisely obtaining specimens of a 
particular type from specific regions of interest of the 
joint have been described previously [15]. The tissues 
were manipulated at all times under liquid nitrogen to 
prevent RNA degradation. Regions of interest were har- 
vested for histological analysis of the articular cartilage, 
and subchondral bone structural parameters were deter- 
mined by micro-computed tomography; regions imme- 
diately adjacent were harvested for RNA isolation and 
microarray analysis. For microarrays, RNA was isolated 
from subchondral bone from four regions: the outer and 
inner lateral tibial (LT) plateau and the inner and central 
medial tibial (MT) plateau (Additional file 1). The gene 
expression intensities changed little within the LT plateau 
(between inner and outer LT regions) and within the MT 
plateau (inner and central MT regions) (Additional file 2); 
we therefore based the analyses in the current study on a 
comparison of gene expression from subchondral bone of 
20 outer LT plateaus and 20 MT plateaus from knee OA 
specimens. These 40 samples were derived from 30 indi- 
vidual knee OA specimens, yielding 10 paired LT and MT 
samples and 10 unpaired LT and MT samples. The 40 OA 
samples chosen for microarray analysis represented a sub- 
set from our previous study [15], and were chosen to have 
LT OARSI scores <8 and MT OARSI scores >18. In 
addition, five LT and five MT paired subchondral bone 
samples from non-OA donors were used as controls for 
the microarray analysis. The characteristics of the samples 
used for microarray analysis are summarized in Table 1 
and Additional file 3. 

The study was approved by the institutional review 
board of Academia Sinica, Tri-Service General Hospital 
and Taipei Medical University Hospital, and written in- 
formed consent was obtained from all of the participants. 

Subchondral bone tissue harvest and RNA isolation 

All tissues were sectioned to generate osteochondral slabs 2 
to 2.5 mm thick from regions of interest for micro- 
computed tomography scanning and histological analysis 
(Figure 1); regions immediately adjacent were isolated for 
subsequent grinding and RNA isolation using a custom- 
made workstation as described previously [15]. Appro- 
ximately 100 mg powdered subchondral bone from the tib- 
ial plateau was obtained for RNA extraction. Five milliliters 
of Trizol (Invitrogen, Carlsbad, CA, USA) were added to 
100 mg (dry weight) ground tissue powder and mixed by 
vortexing until homogeneous. A total of 1 ml chloroform 
was added to the homogenized sample and centrifuged at 
12,000 for 15 minutes at 4°C to achieve phase sepa- 
ration. The clear aqueous phase was transferred to a 
fresh tube for a second phase separation with phenol/ 
chloroform/isoamyl alcohol (25:24:1; Sigma- Aldrich, St. 
Louis, MO, USA), followed by RNA precipitation with pure 



Chou et at. Arthritis Research & Therapy 2013, 15:R190 
http://arthritis-research.conn/content/1 5/6/R1 90 



Page 3 of 1 2 



Table 1 Characteristics of the samples used for microarray analysis 



Osteoarthritis 



Nonosteoarthritis 



Characteristic 


LT 


MT 


LT 


MT 


OARSI 


4.7 ± 1 .34 


Z_w.U J_ 1 .vJZ_ 


J.WW J_ 1 .(J-J 


3.75 ± 1 .67 


BV^ (mm^) 


21.32 ±7.47 


65.39 ±12.61* 


24.41 ±3.71 


35.26 ±10.1 2 


SMI 


1 .56 ± 0.34 


-1.69 ±1.70* 


0.99 ±0.21 


0.68 ± 0.73 


Tb.Th (mm) 


0.1 8 ±0.04 


0.36 ± 0.09* 


0.1 9 ±0.01 


0.26 ± 0.08 


Tb.N (1/mm) 


1.16 ±0.27 


1 .83 ± 0.33* 


1.28 ±0.21 


1 .36 ± 0.09 


Tb.Sp (mm) 


0.51 ±0.12 


0.25 ± 0.09* 


0.56 ± 0.09 


0.49 ± 0.03 


Age (years) 


71.10 ±9.55 


69.65 ± 9.55 


38.40 ±13.45 


38.40 ±13.45 


Male sex (%) 


33 


34 


36 


32 


Number 


20 


20 


5 


5 



Values presented as mean ± standard deviation. BV/TV, percent bone volume; LT, lateral tibial; MT, medial tibial; OARSI, histological scoring system to grade 
severity of osteoarthritis; SMI, structure model index; Tb.N, trabecular number; Tb.Sp, trabecular separation; Tb.Th, trabecular thickness. *P <0.001 considered 
significant comparing MT osteoarthritis versus LT osteoarthritis. P values were determined by two-tailed unpaired Student's f test. 



isopropanol. After washing with 70% ethanol, the air-dried 
pellet was dissolved in RNase-free water and genomic 
DNA was removed by DNase digestion (Qiagen, Hilden, 
Germany) according to the manufacturers instructions. 
To further improve RNA quality, one-half volume of 
7.5 M LiCl (Invitrogen) was added and incubated at -20°C 
for at least 30 minutes followed by RNA precipitation 
by centrifugation at 12,000 for 15 minutes. After 
70% ethanol washing, the RNA was resuspended in 
RNase-free water. The RNA concentration and purity 
were determined by Nano-Drop (NanoDrop Technologies, 
Wilmington, DE, USA). The RNA integrity number 
and the 28 s/18 s ratio were estimated using RNA 6000 



Nano Assays on an Agilent 2100 Bioanalyzer (Agilent 
Technologies, Santa Clara, CA, USA) according to the 
manufacturer s instructions. 



Microarray analysis 

Four hundred nanograms of total RNA per sample were 
used for one round of cRNA synthesis and amplification. 
Cyanine 3-labeled cRNA was then purified and hybri- 
dized to Agilent whole human genome 44 k microarray 
chips (Agilent Technologies). All procedures were car- 
ried out according to the manufacturers instructions 
(Agilent Technologies, Taipei, Taiwan). The array signal 
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Figure 1 Sites of micro-computed tomography analysis and regions chosen for RNA isolation. (A) Human osteoarthritic knee tibial plateau 
(osteoarthritis (OA)) and (B) normal knee tibial plateau (non-OA). Red dashed boxes, regions chosen for RNA isolation; black boxes, regions 
chosen for micro-computed tomography and histological analysis at the two regions of interest. (C) Representative three-dimensional images of 
subchondral bone from the lateral tibial (LT) and medial tibial (MT) plateau from OA and non-OA specimens. 
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intensities were further analyzed by the Agilent Gene- 
Spring GX software (version 11.5; Agilent Technologies). 

Gene expression values of all datasets were normalized 
using quantile normalization; probes with low signal in- 
tensities were excluded by setting the filter above 32. 
The normalized values were log transformed. Significant 
differentially expressed genes between lateral and medial 
samples in OA patients were identified using an adjusted 
Students t test {P <0.05) corrected for multiple compa- 
risons with the Benjamini-Hochberg test. The results 
were used to run principal component analysis to un- 
cover the trends among OA and normal samples. Dif- 
ferentially expressed genes between medial and lateral 
regions among groups were identified by a threshold 
of >2 fold-change and P <0.05. The hierarchical clus- 
tering method, with Euclidean distance and centroid 
linkage, was further applied to cluster the differences in 
gene expression levels between samples. Lists of dif- 
ferentially expressed genes were imported into ingenuity 
pathway analysis (IPA; Ingenuity Systems, Redwood City, 
CA, USA) to identify functional annotations and bio- 
logical interactions from the many known relationships 
between proteins, genes and diseases. The biological 
interaction scores were defined by the IPA statistical 
algorithm based on its knowledge base, and the adjusted 
P value was calculated by Fishers exact test and cor- 
rected for multiple comparisons with the Benjamini- 
Hochberg test. All microarray raw data are available 
through the GEO database [GEO:GSE51588]. 

Quantitative reverse-transcription polymerase chain 
reaction validation 

To validate results from microarray analysis, we per- 
formed reverse-transcription polymerase chain reaction 
(qRT-PCR) for 85 genes (plus GAPDH) on nine additional 
OA paired LT and MT subchondral bone specimens; eight 
of these additional OA knee joint specimens have been 
described previously [15], as well as the qRT-PCR results 
for some of the genes combined with data yielded from 
one additional knee OA specimen. Forty-three of these 
genes were identified for verification on the basis of sig- 
nificant fold-change (>2) comparing LT and MT micro- 
array results. The remaining 42 genes were identified on 
the basis of their potential involvement in OA as des- 
cribed previously [15]. The expression of these 42 genes 
was previously analyzed for seven regions of interest from 
eight OA cartilage and subchondral bone specimens. 

In the current study, we analyzed two regions of interest 
(LT and MT subchondral bone corresponding to outer LT 
and center MT of the prior work) in nine specimens (the 
eight former specimens and an additional OA specimen 
with available LT and MT regions). The qRT-PCR was per- 
formed using the Taqman high-density microfluidic cards 
(Invitrogen, Carlsbad, CA) according to the manufacturer s 



instructions. RNA from each region was converted into 
cDNA using Superscript III reverse transcriptase (Invitro- 
gen, Carlsbad, CA). The qRT-PCR reactions were per- 
formed on the ABI Prism 7900HT Sequence Detection 
system and the fluorescent signal intensity was analyzed 
by Sequence Detector v2.3 software (Applied Biosystems, 
Foster City, California). The levels of cDNA among sam- 
ples were normalized to the expression of GAPDH and an- 
alyzed with the AACt relative quantification method to 
identif)^ significant variation in gene expression comparing 
LT and MT regions in OA samples. Two-tailed unpaired 
Student t tests were performed to evaluate for statistically 
significant differences in gene expression levels between 
regions. 

Immunohistochemistry staining 

Approximately 2 mm diameter sections from the tibial 
plateau were fixed in 4% paraformaldehyde (Sigma- 
Aldrich) overnight and decalcified in 10% ethylenedia- 
mine-A/, A/, A/", A/"-tetraacetic acid, disodium salt, dehydrate 
(Sigma- Aldrich) for 2 weeks. After decalcification, the sec- 
tions were embedded in paraffin and 10 \im sections were 
prepared. The embedded tissue sections were then depa- 
raffmized, hydrated, and the endogenous peroxidase activity 
was quenched in 0.3% H202-methanol for 30 minutes. 
After incubating in 0.5% trypsin/ 1% CaCl2 at 37°C for 
30 minutes and blocking with diluted goat serum at 25°C 
for 1 hour, the sections were incubated with anti-human 
PSOTN (dilution 1:300) antibodies (abl4041; Abeam, 
Cambridge, MA, USA), and LEP (dilution 1:1,000) anti- 
bodies (ab 16627; Abeam) or nonimmune goat serum at 
4°C overnight and reacted with biotinylated secondary anti- 
body for 30 minutes. Signal amplification and staining were 
performed using Vectastain ABC kits (Vector Laboratory, 
Burlingame, CA, USA) according to the manufacturers 
protocol and counterstained with hematoxylin solution Gill 
no. 3 (Sigma- Aldrich). 

Results 

Analysis of differential gene expression 

For the microarray analysis, 972 differentially expressed 
genes with >2 fold-changes between OA MT and OA LT 
regions were identified; 420 of these genes were upregu- 
lated and 552 were downregulated in OA MT compared 
with OA LT (Additional file 4). In total, 640 genes were 
differentially expressed with greater than twofold expres- 
sion level changes (308 upregulated and 332 downre- 
gulated), 163 genes with changes greater than threefold 
(65 upregulated and 98 downregulated). 111 genes with 
changes greater than fourfold (28 upregulated and 83 
downregulated), 35 genes with changes greater than 
sixfold (11 upregulated and 24 downregulated), and 23 
genes with changes greater than eightfold (eight upregu- 
lated and 15 downregulated) (Figure 2A). A graphical 
depiction of microarray gene expression results across 
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Figure 2 Microarray analyses of gene expression of osteoarthritis and non-osteoarthritis of subchondral bone and qRT-PCR validation 
in osteoarthritis knee subchondral bone. (A) Differentially expressed genes were classified based on their expression levels with a minimum of 
twofold, threefold, fourfold, sixfold, and eightfold changes. (B) Unsupervised hierarchical clustering of osteoarthritis (OA) and non-OA samples was 
performed for genes whose differential expression exceeded twofold (972 genes). Distances between samples were detected with a Euclidean 
algorithm and clustered with a centroid linkage method. The OA medial tibial (MT) samples (gray) were clustered together and clearly separated 
from the OA lateral tibial (LT) (brown) and non-OA samples. Non-OA LT samples (N-LT, red) and non-OA MT samples (N-MT, blue) were clustered 
together with OA LT samples. (C) Eighty-five genes (plus GAPDH) were selected for quantitative reverse-transcription polymerase chain reaction 
(qRT-PCR) validation on a separate group of nine OA subchondral bone specimens. The results of qRT-PCR were strongly consistent with those of 
microarray analysis (/?^ = 0.5764, P <0.0001). PC, fold-change. 



the entire extent of the tibial subchondral bone can be 
seen in Additional file 4. 

Unsupervised hierarchical clustering was performed 
on the differentially expressed genes to identify the 



groups or clusters of genes based on similarities between 
the expression profiles of the samples. On the basis of 
gene expression profiles, OA LT and OA MT samples 
were grouped into two distinct clusters. However, the 
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pattern for LT and MT samples from non-OA donors 
was mixed and did not cluster by sample site; the clus- 
tering of all non-OA samples was more akin to the OA 
LT samples than to the OA MT samples (Figure 2B). 

Quantitative reverse-transcription polymerase chain 
reaction validation 

Nine joint tissue samples were available for validation of 
microarray results; qRT-PCR analysis of LT and MT sub- 
chondral bone regions of these samples was performed on 
86 selected genes with GAPDH as endogenous control 
Based on the results of histological evaluation and micro- 
computed tomography scanning, the subchondral bone 
samples selected from the OA MT regions displayed cha- 
racteristics of severe stages of OA; whereas the OA LT 
regions displayed characteristics of relatively normal 
joints. The mean expression pattern of the LT region was 
compared with the mean expression pattern of MT re- 
gions. By qRT-PCR, 77 of 85 genes showed a statistically 
significant difference between regions (Additional file 5). 
Expression by qRT-PCR matched the pattern (upre- 
gulation or downregulation) of expression observed by 
microarray, and the two datasets were strongly correlated 
{R^ = 0.5764, P <0001; Figure 2C). The most differentially 
upregulated genes in OA MT regions were ILll, POSTN, 
ASPK COL6A3, TUBB3, C0L3A1, OGN, DI02, PPEFl, 
and TNFSFll, The most differentially downregulated 
genes in OA MT regions were LEP, APOB, ADHIB, 
CCL8, TUSCS, KLB, NPR3, PLINl, PCKl, and MY016. 
All of these genes were validated by qRT-PCR except for 
ADHIB and MY016y because these two genes were not 
available on the microfluidic cards used for qRT-PCR 
(Table 2). 

Functional and pathway analysis 

IPA of gene expression profiles revealed 279 functional 
networks (canonical pathways) (Additional file 6) that 
were significantly associated with the differentially ex- 
pressed genes. The predominant associated networks and 
biological functions were related to lipid metabolism and 
mineral metabolism (Additional files 7 and 8), connective 
tissue disorders, cellular growth and proliferation, connect- 
ive tissue development and function (Additional file 9). 
Most of these pathways were associated with OA patho- 
genesis such as bone remodeling, adiposity, connective 
tissue, and cell activity and death (Table 3). Interestingly, 
on performing pathway analysis with the most upregu- 
lated and downregulated genes, a novel signaling network 
in OA was identified; this network was generally cha- 
racterized by upregulation of bone mineralization or 
collagen-associated genes, such as C0L3A1 (upregulated 
7.33-fold, ranking sixth in the top upregulated gene list), 
COL6A3 (upregulated 8.27-fold, ranking fourth in the top 
upregulated gene list), BMPl, BMP7, POSTN (upregulated 



15.65-fold, ranking second in the top upregulated gene 
list), WISPh HTRAh SOST, and ITGAll, and the down- 
regulation of genes associated with cellular metabolism, 
proliferation or differentiation {NMB, LEP (downregulated 
23.27-fold, ranking first in the top downregulated gene 
list), CHRDL2, GRB14, CIDEC, CIDEA, BTG2, PLAC8, 
and PRDM16), Additional differentially expressed genes in 
this network included ANGPTLl (upregulated 2.16-fold), 
which is involved in angiogenesis. Both POSTN and LEP 
may affect the phosphorylation of AKT signaling 
(Figure 3A). 

Validation of protein expression in subchondral bone 
specimens 

The regions of interest of subchondral bone samples 
ground for the gene expression study were about 5 mm 
deep under the calcified cartilage, which included the sub- 
chondral bone plate and subarticular spongiosa. These 
regions have a very heterogeneous complement of cells, 
including osteoblasts, osteoclasts, osteocytes and bone 
marrow cells. Although the specific cell types contributing 
to the changes in gene expression cannot easily be con- 
firmed, all of the cell types in the subchondral bone would 
be expected to contribute to the subchondral bone gene 
expression profile. To obtain further independent valid- 
ation of the microarray results, immunohistochemical 
staining was performed to investigate the protein expres- 
sion of proteins Periostin and Leptin corresponding to 
genes POSTN and LEP in human OA and non-OA tissue 
sections from LT and MT regions. The highly expressed 
Periostin could be found in osteoblasts beneath the tide- 
mark, in osteocytes, and in lining cells of the OA MT re- 
gions, but not OA LT regions (Figure 3B,C). Protein 
expression of Leptin was low in OA MT regions, but 
strong in OA LT regions (Figure 3D,E). By immunohisto- 
chemical analysis, the proteins Periostin and Leptin were 
not differentially expressed in the MT and LT regions of 
non-OA samples (Additional file 10). 

Discussion 

To the best of our knowledge, this represents the first 
comprehensive whole genome expression profiling of 
non-OA and OA subchondral bone. The validity of the 
microarray results was confirmed by qRT-PCR for a se- 
lected subset of genes. Based on analysis by IPA, we 
identified biological functions implicated in the patho- 
genesis of OA in bone. Many of the functional networks 
identified in this study were also associated with pro- 
cesses in connective tissue development and function. 

In our previous study [15], we evaluated the expres- 
sion of 61 OA genes and correlated the expression pat- 
terns with the bone morphometric measurement (bone 
volume). The expression of a total of 45 of 61 genes in 
the subchondral bone was significantly correlated with 
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Table 2 Top ten upregulated and downregulated genes comparing lateral tibial bone (below intact cartilage) with 
medial tibial bone (below lesioned cartilage) 



Accession 
number 



Gene 
symbol 



Gene name 



Fold-change (Ml vs. LI) 

Pvalue^ OA Non-OA OA 

(microarray) (microarray) (qRT-PCR) 



Genes with increased expression in medial tibial bone 



NM_ 


_000641 


ILll 


lnterleukin-1 1 


1.32x10" 


-14 


16.02 


1.12 


iio.3r 


NM_ 


_006475 


POSTN 


Periostin, osteoblast specific factor 


1.48x10" 


-10 


15.65 


1.55 


9.68* 


NM_ 


_0 17680 


ASPN 


Asporin 


2.82x10 


-9 


9.04 


2.93 


6.34* 


NM_ 


_057165 


COL6A3 


Collagen, type VI, alpha 3 


9.65x10" 


-10 


8.27 


1.33 


5.ir 


NM_ 


_006086 


TUBB3 


Tubulin, beta 3 


6.98x10" 


-15 


7.36 


1.11 


6.28* 


NM_ 


_000090 


C0L3A1 


Collagen, type III, alpha 1 


5.33x10" 


-16 


7.33 


1.32 


6.45* 


NM_ 


_033014 


OGN 


Osteoglycin 


7.24x10" 


-10 


6.69 


1.70 


9.33* 


NM_ 


_0 13989 


DI02 


Deiodinase, iodothyronine, type II 


1.58x10" 


-1 1 


6.20 


1.31 


5.74* 


NM_ 


_006240 


PPEFl 


Protein phosphatase, EF-hand calcium 
binding domain 1 


5.31 xlO" 


-13 


6.14 


1.97 


4.51* 


NM_ 


_003701 


TNFSFl 1 


Tumor necrosis factor (ligand) superfamily, 
member 1 1 


1.21 xlO" 


-1 1 


5.90 


1.07 


2.64* 



Genes with decreased expression in medial tibial bone 

NM_000230 LEP Leptin 

NM_000384 APOB Apolipoprotein B 

NM_000668 ADHIB Alcohol dehydrogenase 1 B (class I), 

beta polypeptide 

NM_005623 CCL8 Chemokine (C-C motif) ligand 8 

NM_1 72367 TUSC5 Tumor suppressor candidate 5 

NM_1 75737 KLB Klotho beta 

NM_001 204375 NPR3 Natriuretic peptide receptor C 

NM_002666 PLINl Perilipin 1 

NM_002591 PCKl Phosphoenolpyruvate carboxykinase 1 

NM_001 198950 MY016 Myosin XVI 



4.00x10"'° 
3.31 xlO"'° 
2.27x10"^ 

8.17 X 10""^ 
4.22x10"'° 
2.24 X 10"^ 
1.36x10"'° 
2.98x10"^ 
1.03x10"'' 
5.81 X 10"^ 



-23.27 
-16.88 
-12.23 

-10.50 
-10.25 
-10.00 
-9.95 
-9.60 
-9.01 



-2.05 
-2.12 
-1.95 

-1.59 
-1.89 
-1.63 
-1.38 
-1.81 
-1.87 
-1.93 



-18.99* 
-10.20* 
NA 

-6.00* 
-8.19* 
-6.09* 
-7.81* 
-5.88* 
-5.44* 
NA 



LI, lateral tibial; MI, medial tibial; NA, not analyzed; OA, osteoarthritis. ^Adjusted P value for microarray corrected for multiple testing by the Benjamini-Hochberg 
method. <0.001 for quantitative reverse-transcription polymerase chain reaction (qRT-PCR), comparing OA MT versus OA LT. P values were determined by two- 
tailed unpaired Student's f test. 



the alteration of bone structure including ADAMTSl, 
ASPN, BMP6, BMPER, CCL2, CCL8, COLSAl, COL6A3, 
C0L7A1, C0L16A1, FRZB, GDFIO, MMP3, OGN, OMD, 
POSTN, PTGES, TNFSFll, WNTl and others, including 
the ratio of OPG (also known as TNFRSFllB) to RANKL 
(also known as TNFSFl 1), RANKL is primarily produced 
from osteocytes, osteoblasts and/or marrow stromal cells 
and is a key factor for osteoclast differentiation and acti- 
vation. RANKL binds to OPG, also known as osteo- 
clastogenesis inhibitory factor, and interacts with the 
receptor (RANK) for RANKL, to modulate the level of 
osteoclast activity and regulate bone homeostasis in re- 
sponse to various endocrine, paracrine, cytokine and 
mechanical stresses [17,18]. As shown in Additional file 11, 
significant correlations were observed for OPG/RANKL 
and the bone parameters including percent bone volume, 
bone structure (structure model index), trabecular thick- 
ness and trabecular space. These data confirmed our 



hypothesis that investigations of subchondral bone gene 
expression changes could provide clues to the pathogenic 
mechanisms of OA and inform development of new diag- 
nostic markers and therapeutic targets. We focused on the 
outer LT and central MT regions for the purposes of com- 
prehensive microarray analyses in the current study. The 
characteristics of the central MT region generally included 
cartilage matrix loss, cyst formation within cartilage 
matrix, denuded sclerotic bone, and thickened subchon- 
dral bone plate, consistent with a bone sclerosis phenotype 
[19], whereas the outer LT region was relatively normal. 
Many of the most significantly regulated genes identified 
by this microarray study have documented roles in the 
pathogenesis of OA, arthritis or bone formation including 
seven of the top 10 upregulated genes {POSTN, ASPN, 
COL6A3, C0L3A1, OGN, DI02, TNFSFll) and two of 
the top 10 downregulated genes [LEP, APOB) [20-27]. 
Through IPA we could also identif)^ the involvement of a 
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Table 3 Genes representative of gene clusters associated with osteoarthritis 



Functional annotation 
(number of genes, P value)^ 






Representative genes 




Fold-change (MT vs. LT) 




Accession 
number 


Gene 
symbol 


Gene name 


P value^ 


OA 
(microarray) 


Non-OA 
(microarray) 


OA 
(qRT-PCR) 


Abnormal morphology 
of bone 


NM_ 


_004370 


COL 12A 1 


Collagen, type XII, alpha 1 


6.55 X 1 0"°^ 


2.26 


1 .01 


1 .76 


(42 genes, P=].5x 10"^) 


NM_ 


_000594 


TNF 


Tumor necrosis factor 


1.17 X 10"°^ 


-3.17 


-1.45 


-3.42* 


Abnormal morphology 
of collagen fibrils 


NM_ 


_000093 


COLS A 1 


Collagen, type V, alpha 1 


2.63 X 1 0"°^ 


2.90 


1.35 


6.40* 


(7 genes, P= 1.68x10"^) 


NM_ 


_000393 


COL5A2 


Collagen, type V, alpha 2 


2.34x10"^° 


3.12 


1.40 


3.18 


Adiposity (24 genes, P <10"^) 


NM_ 


_003394 


WNTIOB 


Wingless-type MMW 
integration site Family, 

1 1 Icl 1 lUcI 1 UD 


9.76x10"°'' 


2.05 


-1.08 


NA 




NM_ 


_001279 


riDFA 


L-cM (JcaLll IllUUCIliy UrrA IIKc 

effector a 


2.74x10"°^ 


— O.ZD 


1 71 
— 1 ./ 1 


MA 
In A 


Bone mineral density 


NM_ 


_025237 


SOST 


Sclerostin 


z.jy X 1 u 


2.51 


1.44 


4.77* 


[z \ genes, r — i.i x lu j 


NM_ 


_005357 


LIPE 


1 in3<;p l^nrmnnp-<;pn<;iti\/p 


4.40 X 1 U 


-8.71 


-1 .73 


-4.21* 


Angiogenesis 


NM 


004673 


A\/VUr / L / 


A r^/^ i r^/^i o"t"i r^- i I/O 1 
Al lyiUpUlcLII 1 MKc 1 


2.55x10"°^ 


9 1 ^ 
Z. 1 o 


1 n? 
1 .uz 


MA 
In A 


(69 genes, P <10"^) 


NM_ 


_002982 


CCL2 


Chemokine (C-C motif) 
ligand 2 


3.08x10"°^ 


-3.76 


-1.31 


-4.87* 


Osteoarthritis 


NM 


001 851 


C0L9A1 


Collagen, type IX, alpha 1 


1.42x10"°^ 


2.19 


-1.05 


2.35 


(23 genes, P<10"^) 


NM 


001463 


FRZB 


Frizzled-related protein 


1.36x10"^° 


-3.69 


1.38 


-7.72* 


Proliferation of cells 


NM 


001004439 


ITGAll 


integrin, alpha 1 1 


4.05x10"^° 


2.31 


-1.03 


NA 


(303 genes, P <10"^) 


NM 


021 077 


NMB 


Neuromedin B 


1.14x10"°^ 


-3.00 


-1.36 


NA 




NM_ 


_004490 


GRB14 


Growth factor receptor- 

UUUIIU piULclll 1^ 


3.38x10"^^ 


-3.20 


-1.18 


NA 




NM_ 


_006763 


BTG2 


BTG family, member 2 


9.11 xlO"°^ 


-3.67 


-1.31 


NA 


Quantity of connective tissue 


NM_ 


_003474 


ADAM 12 


ADAM metallopeptidase 
domain 12 


1.86x10"^° 


2.30 


1.10 


1.93 


(48 genes, P <10"^) 


NM_ 


_022094 


CIDEC 


Cell death-inducing DFFA-like 
effector c 


2.08x10"°^ 


-6.09 


-1.63 


NA 


Differentiation of cells 


NM_ 


_054034 


FNl 


Fibronectin 1 


1.44x10"°^ 


2.19 


1.38 


2.10* 


(175 genes, P<10-^) 


NM_ 


_016619 


PLAC8 


Placenta-specific 8 


1.01 xlO"°^ 


-2.28 


1.14 


NA 


Cell movement 


NM_ 


_000094 


C0L7A1 


Collagen, type VII, alpha 1 


2.25x10"°'' 


2.28 


-1.00 


3.81* 


(176 genes, P <10"^) 


NM_ 


_012193 


FZD4 


Frizzled family receptor 4 


3.62x10"^° 


-3.02 


-1.37 


4.93* 


Necrosis 


NM_ 


_080838 


WISPl 


WNTl inducible signaling 
pathway protein 1 


3.79x10"°^ 


2.50 


1.09 


NA 


(218 genes, P<10"^) 


NM_ 


_004797 


ADIPOQ 


Adiponectin, CIQ and 
collagen domain containing 


1.71 xlO"°^ 


-6.63 


-1.69 


10.51* 


Mineralization of cells 


NM_ 


_001719 


BMP7 


Bone morphogenetic 
protein 7 


1.65x10"°'' 


2.50 


1.34 


2.26** 


(14 genes, P = 3x10"^) 


NM_ 


_0 15424 


CHRDL2 


Chordin-like 2 


1.40x10"°^ 


-4.04 


1.38 


-15.71* 



LT, lateral tibial; MT, medial tibial; NA, not analyzed; OA, osteoarthritis. ^Number of genes associated with the particular ingenuity pathway and differentially 
expressed in MT versus LT regions of OA samples by microarray analysis; P value calculated by Fisher's exact test. "^Adjusted P value for microarray corrected for 
multiple testing by the Benjamini-Hochberg method. *P <0.001 and **P <0.05 for quantitative reverse-transcription polymerase chain reaction (qRT-PCR); P <0.05 
was considered significant comparing OA MT versus OA LT. P values determined by two-tailed unpaired Student's f test. 



lipid metabolism protein network and a mineral metabo- 
lism protein network that involved 29 and 22 differentially 
expressed genes respectively; these networks included 
PLINE, LIPE, DGAT2, ADRBl, NPYIR, HCAR3 and 
P2RRY14 that changed -9.6-fold, -8.7-fold, -7.8-fold, 
-6.2-fold, -7.4-fold, -5.1-fold, and -5.1-fold respectively 



comparing OA MT with OA LT (Additional files 7 and 8). 
These functional networks indicated that the bone cells in 
the OA MT regions were characterized by a condition of 
low energy production that included a low rate of mineral 
metabolism. These results are compatible with a low rate 
of bone remodeling in MT regions. 
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Figure 3 Prediction of an interactive networic of genes and validation of protein expression in osteoarthritis subchondral bone. 

(A) Network analysis was performed on differentially expressed genes in the osteoarthritis (OA) samples by ingenuity pathway analysis (IPA). 
Relationship labels: A, activated; P, phosphorylation/dephosphorylation. (B), (C), (D), (E) Antibodies against proteins Periostin and Leptin 
corresponding to genes POSTN and LEP were used to detect protein expression (dark brown staining) by immunohistochemistry of OA 
subchondral bone. Highly expressed, Periostin could be found in osteoblasts (OB) under the tidemark (TM), in osteocytes (OC) and in lining cells 
(LC) in the OA medial tibial (MT) regions, but not in the OA lateral tibial (LT) regions. Leptin could barely be detected in the OA MT regions, but 
was strongly expressed in the OA LT regions. All sections were counterstained light blue with hematoxylin. Scale bar = 0.2 mm. 



Increased bone remodeling is an important predictor of 
OA progression [8] and is characteristic of early stages of 
OA including resorption of bone and increased porosity in 
the subchondral bone region [28]. Given that OA MT re- 
gions of interest could be defined as representing late stages 
of OA, we were interested to determine whether OA LT re- 
gions of interest could be described as representing early 
stages of OA. The OARSI histopathological scores of all OA 
LT samples were <6, representing a relatively intact overlying 
cartilage matrix with minimal superficial fibrillation not too 
dissimilar from normal. However, in contrast to normal sam- 
ples, the bone structural parameters of the OA LT regions 
had a lower bone volume, thinning of the subchondral bone 
plate, increased porosity in the subchondral region and re- 
duced bone density (bone density data not shown) suggestive 
of early OA. Another indicator of early OA, as demonstrated 
by animal models of OA and in OA patients, is increased 
bone remodeling characterized by a decreased expression ra- 
tio of OPG to RANKL [29,30]. In contrast to non-OA tissues, 
a reduced OPG to RANKL ratio was identified in OA LT re- 
gions compared with non-OA LT regions {OPG lacked sig- 
nificant change, but RANKL was upregulated 1.63-fold, one- 



tailed unpaired Student t test P = 0.048), consistent with 
early OA. Thus, based on the bone parameters and OPG 
to RANKL ratios, comparing LT with non-OA we con- 
cluded that the bone gene expression patterns of OA LT 
regions were consistent with an early-stage OA phenotype. 

According to a hypothetical model of OA pathogenesis 
proposed by Burr and Gallant [6], joint loading increases 
bone remodeling in the early stage of OA, and a bone- 
cartilage crosstalk may occur via diffusion of small mole- 
cules or increased vascularization at the deep layers of 
cartilage that could interfere with the normal collagen net- 
work. These changes would lead to a loss of aggrecan, in- 
creasing inflammation, and result in bone sclerosis at late 
stages of OA. Several of the genes identified in this study 
fit nicely into this hypothetical model. For instance, upre- 
gulation of SOST was identified (2.51-fold); SOST is an 
inhibitor of Wnt signaling that responds to mechanical 
loading in the proximal rat tibia and is associated with 
bone mass changes [31]. In addition, several genes rele- 
vant to cell differentiation and activity of osteoblasts, 
osteoclasts and osteocytes were identified in both a mouse 
load-regulated gene expression model [32] and in this 
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human study; examples included ASPN (regulator of 
osteoblast collagen mineralization), WNT5A (agonist of 
WNT signaling pathway in osteoblasts and osteocytes) 
and VCAN (regulating transforming growth factor beta 
expression in osteoclasts). These findings suggest that 
mechanical loading plays an important role in structural 
changes of human subchondral bone, initiates bone re- 
modeling, and contributes to the pathogenesis of OA. 
Moreover, our results could provide a possible novel mo- 
lecular mechanism to explain changes in bone remodeling 
during OA development based on an AKT-regulated net- 
work. AKT is a serine/threonine-specific protein kinase 
that plays a pivotal role in many cellular processes such as 
metabolism, apoptosis, cell proliferation, transcription, cell 
migration and extracellular matrix alternation [33]. In an 
early-stage of OA, overexpressed LEP may increase phos- 
phorylation of AKT [34], which will trigger the down- 
stream signaling pathway to increase bone remodeling. In 
a late stage of OA, upregulated POSTN may inhibit phos- 
phorylation of AKT [35] and will decrease the cellular me- 
tabolism, cell proliferation, and differentiation; the rate of 
bone remodeling would thereby decline. 

One of the limitations of this study was that the sub- 
chondral bone samples were obtained from OA joints at 
end-stage disease due to the difficulty of obtaining joints 
with early stages of OA from humans. However, compa- 
ring our results with the gene expression profiles of sub- 
chondral bone in an early OA experimental mouse model 
[12], their reported fold-changes in SB gene expression of 
ASPN, CCL2, COL3A1, COL5A1, POSTN and TNFSFll 
at the initial stages of OA were similar to those of our 
study. This result strongly supports the credibility of this 
study. Another limitation was the difficulty distinguishing 
whether differentially expressed genes reflected OA pro- 
gression in cartilage, differences in innate bone structure, 
or were driven by changes in mechanical loading. In our 
currently study, significant correlation between overlying 
cartilage integrity and underlying subchondral bone struc- 
tural parameters could be identified (Additional file 12), 
suggesting that differential gene expression across the sub- 
chondral bone is due, at least in part, to differences in 
cartilage integrity of the LT and MT regions in OA. 
Microarray analysis of additional non-OA control samples 
could help to exclude differentially expressed genes driven 
by innate bone structure. However, to exclude genes 
driven by mechanical loading, microarray analysis of 
lateral compartment dominant OA is necessary. 

The site of initiation of OA, at bone or cartilage, has 
been disputed for decades. Numerous clinical and ex- 
perimental studies have confirmed that increased bone 
volume and changes in bone quality of the tibial subchon- 
dral bone of the knee are related to loss of cartilage integ- 
rity [7,8,36-38]. Goldring and Goldring have pointed out 
that subchondral bone responds more rapidly to adverse 



loading and events than cartilage [5]. Moreover, the integ- 
rity of subchondral bone predicts later worsening of radio- 
graphic osteoarthritis progression [8]. As summarized by 
Martel-Pelletier and Pelletier, an abnormal rate of bone re- 
modeling, incomplete bone mineralization and increased 
osteoid collagen matrix will result in hypomineralization of 
subchondral bone and deterioration of cartilage as mani- 
fested by knee OA progression [39]. These data would sug- 
gest that drugs designed to target both the bone and 
cartilage compartments to optimize bone remodeling rates 
and inhibit cartilage loss in the early stages of OA could in- 
hibit progression of disease. The recent success of strontium 
ranelate as a disease-modifying agent for OA provides the 
first clear evidence that this hypothesis has been confirmed 
[40]. Some other bone-targeting agents have been assessed 
as OA therapeutics, including bisphosphonates [41,42], calci- 
tonin [43] and vitamin D [44]. These bone-acting agents are 
unlikely to be effective in late stages of OA, when remodeling 
is already suppressed. The development of diagnostic bio- 
markers to detect progressive changes in early stages of OA 
is therefore important. The genes highly regulated in early 
stages of OA in our microarray data could provide potential 
diagnostic biomarkers or therapeutic targets for OA. 

Conclusion 

Our microarray results were obtained by examining early 
and late stages of OA and non-OA samples. Most of the 
differentially expressed genes in bone are involved in car- 
tilage and bone development, OA pathogenesis and bone 
remodeling in the early and late stages of OA. This direct 
assessment of gene expression profiUng in OA subchon- 
dral bone provides a comprehensive understanding of dis- 
ease pathology that could further the development of new 
OA biomarkers and therapeutic strategies. 

Additional files 



Additional file 1: Shows sites of analysis and regions for RNA 
isolation. 

Additional file 2: Shows self-organizing maps analysis (SOM) of the 
differentially expressed genes across the entire tibial plateau. 

Additional file 3: Presents characteristics of the paired and 
unpaired OA samples used for microarray analysis. 

Additional file 4: Presents a list of differentially expressed genes. 

Additional file 5: Presents gene validation by qRT-PCR. 

Additional file 6: Presents a list of canonical pathways. 

Additional file 7: Shows a potential network of lipid metabolism in OA. 

Additional file 8: Shows a potential network of mineral metabolism in OA. 

Additional file 9: Presents a list of dominant functional networks. 

Additional file 10: Shows validation of protein expression in 
human non-OA knee samples. 

Additional file 11: Presents correlation of OPG/RANKL and 
structural parameters of the subchondral bone. 

Additional file 12: Presents correlations between cartilage integrity 
and the bone parameters. 
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